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' ABSTRACT 

o : 

+3 • We investigate the conditions for the presence of a magnetically inactive dead 

zone in protostellar disks, using 3-D shearing-box MHD calculations including 
. vertical stratification, Ohmic resistivity and time-dependent ionization chemistry 

^ . Activity driven by the magnetorotational instability fills the whole thickness of 

the disk at 5 AU, provided cosmic ray ionization is present, small grains are ab- 
sent and the gas-phase metal abundance is sufficiently high. At 1 AU the larger 
column density of 1700 g cm -2 means the midplane is shielded from ionizing 
particles and remains magnetorotationally stable even under the most favorable 
conditions considered. Nevertheless the dead zone is effectively eliminated. Tur- 
bulence mixes free charges into the interior as they recombine, leading to a slight 
coupling of the midplane gas to the magnetic fields. Weak, large-scale radial 
fields diffuse to the midplane where they are sheared out to produce stronger 
azimuthal fields. The resulting midplane accretion stresses are just a few times 
less than in the surface layers on average. 



Subject headings: circumstellar matter — solar system: formation — stars: for- 
mation — instabilities — MHD 
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INTRODUCTION 



Turbulence plays a central role in protostellar disks. It tra nsports orbital angular mo - 
mentum, controlling the accretion of gas on the central star ( Shakura fc Sunvaevl 1973 ): 



bring s dust grains together in the early stages of planet fo rmation ( Weidenschil 



19931 ); and modifies the orbital migration of protoplanets (INelson fc Papaloizoul 12004 ) . Tur- 



ing &: Cuzzi 



bulence can be driven by the magneto- rotational instability (MRI; Balbus & Hawley 1991) 
in regions where the ionization is sufficient to couple the gas to the magnetic fields. The 
MRI operates by transferring orbital angular momentum outward along magnetic field lines 
linking gas on different orbits. The fastest linear growth rate is three-quarters the orbital 
frequency Q and occu rs at wavelengths near X r = 2ttva/Q in fully-ionized Keplerian disks 
with Alfven speed va ( IBalbus &: Hawleylll99ll ). In the non-linear regime, the gravitational 
energy released from the insp iraling gas sustains turbulence and regenerates the magnetic 
fields ( IBalbus fc Hawley! Il998l ). Good coupling of the fields to the gas within a few tenths 
of an AU of the star co mes from thermal ionization of the low-ionizat ion-potential elements 
potassium and sodium (IPneuman fc Mitchelj|l965l ; IUmebayashilll983l ). At greater distances, 
the ionization is due to the absorption of stellar X-rays or inters tellar cosmic rays in the disk 
material. The X-rays penetrate to columns of about 10 g cm -2 ( llgea fc Glassgoldlll999l ) and 
the cosmic rays to 100 g cm -2 , while recombination is generally fastest in the dense interior, 
leading to a layered structure. The top and bottom surf aces of the dis k are magnetically 
active. The interior is quiescent and called the dead zone (lGammielll996l ). The dead zone is 
of special interest because it overlaps the region where the planets are found and potentially 
provides the environment and the raw ingredients for planet formation. The simplest lay- 
ered picture, however, yields an accretion rate independent of stellar mass, conflicting with 
data for T Tauri stars. Including the extra heat ing from irradiation o f the disk by the star 
gives a mass dependence that is still too weak (IHartmann et al.ll2006l ). indicating that our 
understanding of the dead zone is incomplete. 

The size of the dead zon e depends on the disk surface density , the abundance of the dust 
that hasten s recombination (IWardle fc Nall999i ; ISano et al.ll2000l ). the presence of long-lived 



metal ions (IFromang et a 



2002 



]) and the rate at which the turbulence mixes the charged 



species into the interior (jlnutsuka &: Sand 120051 ; lllgner fc Nelson! l2006bl ) . Well-mixed 0.1- 



(Ltm grains with the interstellar dust-to-gas mass ratio in the minimum-mass Solar nebula 
prod uce a dead zone extending to about 20 AU (ISano et al.l |2000| ; iMatsumura fc Pudritz 
20061 ). Even if small grains are absent, the annulus immediately outside the thermally- 
ionized region is quiescent for a wide range of disk models because the high surface densities 
mean cosmic rays and X-rays are absorbed before reaching the midplane. In this paper 
we investigate whether the dead zone can be eradicated under favorable conditions leaving 
the entire disk active. Because cosmic rays penetrate more deeply, their ionizing effects 
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are included despite possible screening by the protostellar wind. Slow recombination is 
ensured by assuming the dust surface area per unit gas mass has declined through settling, 
coagulation and orbital evolution sufficiently that recombination on grain surfaces can be 
neglected. Since the abundance of gas-phase metal atoms is not well known, it is treated 
as a parameter. We focus on conditions near the present or bits of the Earth and Jupiter 
and use the minimum-mass So l ar neb ula (IHayashi et al.lll985l ) as the disk model. We extend 
results from iFleming fc Stone! (120031 ) by making stratified shearing-box MHD calculations 
including time-dependent ionization chemistry. 

The paper has six sections. The ionization-recombination reaction network, disk model 
and MHD methods are described in section [21 The network is used in section [3] to estimate 
the size of the dead zone in local chemical equilibrium. Expansion or contraction of the 
dead zone through turbulent mixing is explored using MHD calculations in section H] while 
section [5] deals with simple expressions for predicting the boundaries of the dead zone. A 
summary and conclusions are in section [6j 



METHODS 



2.1. Reaction Network 



The ionization state of the disk gas is calculated using the simplified reaction network 
developed for molecular cloud studies by lOppenheimer fc Dalgarnol (119741 ) . Five species are 
involved: molecular hydrogen and its ion; magnesium, a representative metal, and its ion; 
and free electrons. These take part in the four reactions 



H 2 

H+ + e- 
H+ + Mg 
Mg + + e" 



H+ + e- 
H 2 

H 2 + Mg H 
Mg + hu. 



(1) 
(2) 
(3) 
(4) 



The ionization (eq. [T]) is balanced by both dissociative recombination of the molecular ion 
(eq. [2]) and molecule-metal charge exchange (eq. [3]) followed by the radiative recombination 
of the metal ion (eq. H|). The rate coefficients are 3 x 10~ 6 /y/T cm 3 s" 1 for dissociative 
recombination, 3 x 10~ 9 cm 3 s _1 for charge exch ange and 3 x 1Q~ U / \/T cm 3 s _1 for radiative 
recombination at a temperature T in Kelvin (jllgner fc Nelson! l2006al ). Because radiative 
recombination is slowest, the positive charge typically resides on metal ions in chemical 
equilibrium. Comparison with a comprehensive network of almost 2000 reactions shows that 
the Oppenheimer & Dalgarno network yields slightly higher equilibrium electron abundances 
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because fewer recombination pathways are included (jllgner fc Nelson! l2006al ). A similar 
reduced reaction set adequately reproduces the ionizatio n fraction from a det ailed chemical 
model of the midplane of a protostellar disk at 1 AU (jSemenov et al.l 120041 ) . The kinetic 
equations for t he reactions [T0 form a stiff set of ODEs that is integrated using semi- implicit 
extrapolation (IPress et al.lll992l ). Solving the reaction network uses about ten times as many 
CPU cycles per zone as evolving the magnetic fields. 

The cosmic ray ionization rate is reduced from the interstellar value (cr — 10~ 17 s _1 



( iCaselli et al.l Il998l ) by the shielding effects of the surrounding disk material. The total 
ionization rate 



C(*) 



CR 



exp 



XCR 



+ exp 



(5) 



XCR 

depends on the cosmic ray absorption depth Xcr — 96 g cm -2 and the column above and 
below the point of interest, x + = f^°° pdz and x~ — fl pdz. In the 3-D MHD calculations 
of section HI the column densities are averaged over the horizontal extent of the domain when 
finding the ionization rate. 



The resistivity is computed from the electron abundance by 



T 

r) = 234 — cmV 1 , 



(6) 



where x e = n e / n n is the ele ctron fraction, n f the electron number density and n n the total 
number density of neutrals (IBlaes fc BalbusHl994Tj . 



2.2. Disk Model 



Our model disk is the minimum-mass Solar nebula in which the surface density and 
temperature vary with radius r as S = 1700(r/AU)~ 3 ' 2 g cm" - 2 and T = 280(r/AU)- 1 / 2 K, 
respectively. The temperature is independent of the distance from the midplane, so that 
in hydrostatic equilibrium the density falls off with height z as exp(— z 2 /2H 2 ). The scale 
height H = c s /fl depends on the sound speed c s = (lZ gas T / (j,) 1 / 2 . 



The composition of the material is Solar (jAsplund et al.ll2005l ). with 0.085 helium and 
3.39 x 10 -5 magnesium atoms per hydrogen nucleus. The hydrogen is molecular and the mean 
molar weight fi = 2.29. Only gas-phase magnesium contributes to the overall ionization level, 
while the magnesium atoms can take three forms, (1) incorporated into silicate and metallic 
grains, (2) adsorbed on grain surfaces and (3) in the gas phase. Me asurements of the diffuse 
interstellar gas show t hat about half the magnesium is in grains ( ISofia. Cardelli fc Savage 
19941 ; lFitzpatrickill997l ) . We assume the remainder is divided between grain surfaces and the 
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gas phase through thermal ad sorption and des o rption and we calculate the equilibrium gas- 
phase fraction as outlined by lllgner fe Nelson! (j2006al ). The equilibrium is reached quickly, 
as the lifetime of adsorbed magnesium atoms is just 10 -4 seconds at 1 AU and 280 K, or 
10 6 seconds at 5 AU and 125 K. Almost all the magnesium not in grain interiors is in the 
gas phase at 1 AU, while at 5 AU just 1.4% stays in the gas if 0.1-micron grains are present 
with the interstellar mass fraction. If the gas-phase magnesium abundance is unchanged by 
the removal of the small grains, its most probable value is half Solar at 1 AU and 0.7% Solar 
at 5 AU. We consider a range of possible magnesium abundances below. 

The evolution of the magnetic fields is potentially affected by the Ohmic, Hall and am- 
bipolar diffusi on terms in the induct ion equation. We estimate the magnitudes of the terms 
as outlined in ISano fc Stond (j2002al ). The Ohmic and Hall terms are comparable near the 
midplane , while the ambipolar term is the large s t of the three only above about 3H. Linear 
analyses (jWardldll999l ; lBalbus fc Terquemll200ll ; ISalmeron fc Wardldl2005l ) indicate the Hall 
term has little effect on the maximum MRI growth rate but extends the range of unstable 
wavelengths. Numerical results suggest the strength of the saturate d turbulence is largel y 



unaffected by the Hall term and is determined by the Ohmic term flSano &: Stone 



2002 



b. 



The a mbipolar term decouples the neutral gas from the ionized component (IHawley fc Stone 
19981 ) at the low densities found in our model disk outside \z\ = 4H, where collisions transfer 
momentum from the ions to the neutrals at rates less than the orbital frequency. Little 
accretion of neutral gas due to MRI turbulence will occur beyond four scale heights. Since 
the Ohmic term is most important for the size of the dead zone, we choose to focus on the 
effects of the Ohmic resistivity. 



2.3. MHD Equations 

The standard equations of resistive MHD, 

^ + V-( P v) = 0, (7) 

_ + v ■ Vv = + (V xB)xB-2fixv + 3fi 2 x - fi 2 z, (8) 

dt p 4vrp v ' y ' 

and 

- = Vx(vxB-i|VxB), (9) 



are sol ved with an isothermal equation of state p = pc 2 s) using the ZEUS code (IStone fc Norman 
1992al lbh. When the chemical reactions are included, an additional continuity equat 



tion 



^ + V-( A v) = 0, (10) 
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is sol ved for each of the five chemical sp ecies i. The domain is a three-dimensional shearing- 
box (IHawley. Gammie. &: Balbud Il995l ) and vertical stratification is included. The radial 
or x-boundaries are periodic, with an azimuthal offset that increases in time according to 
the radial gradient in orbital speed. The azimuthal or y-boundaries are strictly periodic, 
while the vertical or z-boundaries allow outflow but no inflow. The domain extends two 
density scale heights H along the radial direction, eight along the orbital direction and four 
either side of the midplane and is divided into 32 x 64 x f 28 grid zones. Short timesteps are 
avoided by applying a density floor of 1/15 the minimum initial density. The floor typically 
operates in just a few grid zones near domain top and bottom, and for all the calculations 
described in this paper the mass added is much less than the mass escaping through the 
boundaries, which in turn is in most cases less than 1% of the total initial mass per 30 orbits. 
The starting magnetic field has strength B chosen so that the midplane MRI wavelength 
is 0.5H. The components of the field are (B x , B y , B z ) = [0, B sxn(7rx/H), B cos(nx/H)]. 
Magnetic pressure and tension forces initially vanish, as the field strength is uniform and 
the field lines are straight. The net magnetic flux is zero because the field lines rotate 360° 
round the y — z plane across the radial extent of the box. Due to the periodic conditions 
on the side boundaries, the vertical component of the magnetic flux is conserved during the 
calculations. The outflow conditions on the top and b ottom boundaries mean that, unlike 

(120031 ). the radial and azimuthal magnetic 



the calculations discussed by Fleming fc Stone 



fluxes are free to evolve over time. 



3. TIMESCALES 

To gain an overall grasp of the problem, we estimate three timescales affecting the size 
of the dead zone. These are (1) the Ohmic dissipation time for the magnetic fields, (2) the 
turbulent mixing time and (3) the time for ionized gas to recombine and reach chemical 
equilibrium. The dead zone is the region where the magnetic fields dissipate too fast for 
the MRI to grow. The dissipation time is calculated in this section assuming local chemical 
equilibrium so that the electron abundance at each height results solely from local ionization 
and recombination processes. However the dead zone boundaries will shift if ionized gas 
from the surface layers is mixed to the interior before recombining. The three corresponding 
timescales are calculated as follows. 

1. Like any diffusion timescale, the dissipation and mixing times are ratios of a squared 
length to a diffusion coefficient. The diffusion coefficient r\ for resistive dissipation is 
given by eq. O For the length we choose the scale height H, yielding an upper bound 
on the dissipation time in the disk interior, where the magnetic fields can alternate in 
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sign over distances shorter than the scale height. MRI turbulence is expected when the 
dissipation time is longer than the growth time of the instability. For example if the 
characteristic unstable wavelength A c is one-tenth of the scale height then turbulence 
requires a dissipation time, H 2 /rj, longer than about (0.1)~ 2 = 100 orbits. The diffusion 
time in orbits can be converted to a magnetic Reynolds number (? s /{rfl) by multiplying 
by 2tt. 

2. The mixing time depends on the turbulent diffusion coefficient, whi ch is the product of 



the m ean squared velocity fluctuation and the eddy correlation time (IFromang fc Papaloizou 



20061 ). For a preliminary estimate of the mixing time, these quantities are taken from 



ideal-MHD shearing-box calculations, with the Ohmic resistivity switched off in eq. M 

3. The recombination time is measured in the hydrostatic model disk starting with the 
magnesium and electron abundances equal, the magnesium fully-ionized and the molec- 
ular hydrogen neutral. The reaction network is then integrated until the electron 
abundance reaches equilibrium. 

We consider conditions at 1 and 5 AU in the minimum-mass Solar nebula. The mixing times 
are found using turbulent diffusion coefficients from the ideal-MHD calculations Ml and M5 
listed in table [H The coefficients are averaged over the horizontal extent of the domain and 
over a period from 10 to 100 orbits after the turbulence is well-established. 



3.1. Timescales at 1 AU 

Results for 1 AU are shown in figure [TJ Because we aim to find the size of the active 
region under the most favorable conditions, we assume here that the gas-phase magnesium 
abundance is Solar. Nevertheless the magnetic fields dissipate in just three orbits at the 
midplane if the ionization is in equilibrium. The M RI is absent near z = as all modes 



fitting in the disk thickness (jSano fc Miyamal Il999l ) are damped. The recombination time 



is several hundred years at the midplane and also at four scale heights, but is only thirty 
years at two scale heights, where the absolute number density of free electrons is greatest. 
At all heights, the turbulent mixing from the ideal-MHD calculation is faster than the re- 
combination, suggesting mixing could affect the size of the dead zone. The recombination 
and dissipation times depend on the metal abundance as shown in figure [2J The midplane 
Ohmic dissipation time is 3 orbits or less across the entire range of abundances considered, 
indicating that changes in the metal abundance alone will not remove the dead zone. In 
contrast, at z — 2H the dissipation time is longer than 100 orbits and magnetic activity is 
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Fig. 1. — Timescales for Ohmic dissipation (dotted), recombination (dashed) and turbulent 
mixing (solid) as functions of height at 1 AU, assuming Solar abundance for magnesium in 
the gas phase. The mixing timescale is taken from the ideal-MHD calculation Ml. 



expected if the magnesium abundance is greater than 10 -5 Solar. The dissipation and re- 
combination times at all heights are abundance-independent above 10~ 4 Solar, so the results 
with Solar magnesium accurately reflect conditions throughout the likely abundance range: 
in local chemical equilibrium, the surface layers are magnetically active while the midplane 
is dead. 



3.2. Timescales at 5 AU 

Figure [3] shows the height dependence of the timescales at 5 AU. The assumed gas-phase 
magnesium abundance is 1% of Solar, near the equilibrium value for thermal adsorption and 
desorption (section [2T2T) . MRI turbulence is expected at all heights because the dissipation is 
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Fig. 2. — Midplane charged particle abundances (top panel) and timescales (bottom) versus 
gas-phase magnesium abundance at 1 AU. Above, is shown by a long-dash-dotted line, 
Mg + by a short-dash-dotted line and electrons by small circles. Below, the recombination 
timescale is shown by a dashed line and the Ohmic dissipation time per scale height by a 
dotted line. Circles near the left-hand edge mark the values when magnesium is absent. 
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Fig. 3. — Timescales for Ohmic dissipation (dotted), recombination (dashed) and turbulent 
mixing (solid), as functions of height at 5 AU, with a gas phase magnesium abundance 1CT 2 
Solar. The mixing timescale is taken from the ideal-MHD calculation M5. 

slow. The relatively high midplane electron fraction x e = 6 x 10~ 10 is due to faster ionization 
than at 1 AU: cosmic rays reaching z = pass through a column of just 76 g cm~ 2 . The 
turbulent mixing is slower than the recombination in the interior, so the ionization profile 
will be unaffected by mixing. Midplane magnesium abundances greater than 10~ 6 Solar yield 
dissipation slow enough to permit the MRI to develop (figure H]). 



According to the timescale estimates in section [3J layered accretion is expected at 1 AU, 
but turbulent mixing could be fast enough to reduce or dispense with the midplane dead 
layer. At 5 AU, the whole thickness of the disk will be active. Here we explore the magnetic 
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Fig. 4. — Midplane charged particle abundances (top panel) and timescales (bottom) versus 
gas-phase magnesium abundance at 5 AU. Symbols are as in figure EJ above, H^" is shown 
by a long-dash-dotted line, Mg + by a short-dash-dotted line and electrons by small circles. 
Below, the recombination timescale is shown by a dashed line and the Ohmic dissipation 
time per scale height by a dotted line. Circles near the left-hand edge mark the values when 
magnesium is absent. 
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activity using resistive MHD calculations with and without time-dependent ionization. 



4.1. Turbulent Mixing at 1 AU 

In calculation Fl the resistivity is a fixed function of height taken from the equilibrium 
ionization state in the initial hydrostatic disk model. The resulting flow has active surface 
layers and a mi dp lane dead zone (figure [5]). The stresses in the active layers are due mostly to 
magnetic forces, with the contribution from gas pressure gradients three to six times smaller, 
while at the midplane the total stress is an order of magnitude weaker and the magnetic and 
hydrodynamic parts are almost equal. The midplane magnetic stresses result from net radial 
and azimuthal fluxes that build up as fields o f the opposing signs esca pe through the open 



vertical boundaries. Similar calculations by iFleming fc Stone! ( 120031 ) showed still weaker 



midplane magnetic stresses because the vertical boun daries were period i c or i mpermeable 



and the net radial magnetic flux remained zero. Like IFleming fc Stond (120031 ) we observe 
hydrodynamic stresses in the dead zone due to sound waves propagating away from the 
active layers. In run VI the reaction network is solved together with the MHD equations. 
Electrons and ions are carried from place to place in the flow and the resistivity varies with 
all three space coordinates and time. For the first 60 years the flow again consists of active 
surface layers and a midplane dead zone (figures and [6]) • The magnetic fields are weaker 
overall than the ideal-MHD case Ml, leading to slower mixing. Recombination is faster than 
mixing within two scale heights of the midplane (figure [7]), so the resistivity approaches 
its equilibrium profile from figure [1] even when the time-dependent chemistry is included. 
The dissipation time differs most from figure [1] in the active layers, where mixing makes the 
electron fraction and resistivity more nearly uniform. The turbulent mixing has little effect 
on the size of the midplane dead zone. 

The situation changes around 60 years, when a period of stronger turbulence leads to 
faster mixing. Ionization at the midplane increases and the active regions expand toward 
z = 0. The midplane electron fraction peaks at 87 years at x e = 1.2 x 10~ n and then varies 
around 10 -11 till the run ends at 150 years. The flow settles into a new state in which the 
mixing at the midplane is faster than the recombination, but still slower than the dissipation 
of magnetic fields (figure [7]). The fast-growing MRI modes are adequately spatially resolved 
throughout the volume. Outside two scale heights the characteristic wavelength averages 
ten or more grid zones in the vertical direction due to the large Alfven speeds, while nearer 
the midplane the linear modes not damped by the resistivity have vertical wavelengths six 
zones or greater. The mean accretion stresses are similar to the ideal-MHD calculation Ml 
in the surface layers outside \z\ = H, but a few times weaker at the midplane (figure [5]). The 
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Fig. 5. — Mean accretion stress versus height in MHD calculations centered at 1 AU. The 
magnetic stress —B x B y /4:H is shown by solid lines, the hydrodynamic stress pv x Sv y by dotted 
lines; Sv y is the azimuthal velocity relative to the background shear flow. The four panels 
show results from an ideal MHD run (top left), a resistive case with a fixed resistivity profile 
i](z) (top right), and early and late stages from a run with time- varying resistivity (bottom 
left and right, respectively). The early results are averaged from 10 to 60 years, the late 
results from 100 to 150 years. The gas-phase magnesium abundance is Solar. 
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Table 1. Calculation Parameters 



Label 


Physics 


R/AX1 


XMg/Solar 


Duration/orbits 


Ml 


Ideal MHD 


1 




100 


Fl 


Fixed rj(z) 


1 


1 


150 


VI 


Varying r}(x,y,z,t) 


1 


1 


150 


M5 


Ideal MHD 


5 




100 


F52 


Fixed i](z) 


5 


io- 2 


100 


F56 


Fixed r](z) 


5 


1(T 6 


100 


F58 


Fixed r)(z) 


5 


1(T 8 


100 




Time / years 

Fig. 6. — Horizontally-averaged magnetic accretion stress versus height and time at 1 AU 
in the ideal MHD calculation Ml (top panel) and the resistive MHD version with time- 
dependent chemistry VI (two lower panels). The pair of black curves nearest the midplane 
in the middle panel shows where the magnetic Reynolds number Re = c 2 s /{r]VL) = 100, the 
white lines Re = 1000 and the outer black curves Re = 10 4 . Black curves in the bottom 
panel mark where the Lundquist number Lu = v\ z /(r]Q) is unity. The common color scale 
is double logarithmic, with reds for positive stresses, blues for negative stresses and green 
for stresses within 10~ 4 dyn cm -2 of zero. 
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Fig. 7. — Timescales for Ohmic dissipation (dotted), recombination (dashed) and turbulent 
mixing (solid) versus height at 1 AU in the MHD calculation VI with time-dependent ion- 
ization. The results are averaged from 10 to 60 years (top panel) and from 100 to 150 years 
(bottom panel). 



midplane is magnetized and has a significant total accretion stress that is 0.2% of the local 
gas pressure. However the magnetic activity differs in character from the surface layers and 
is also unlike the midplane in the ideal-MHD run Ml in several respects: 



1. The magnetic field shows few features smaller than a quarter of a scale height, as 
expected from the dissipation timescale of 20 years per scale height. 



2. The magnetic stress m xy = —B x B y /4:7i is due mostly to large-scale fields, similar to 
the dead zone during the early part of the calculation, rather than small-scale fields 
as in the ideal-MHD run (figure [8]). The horizontally-averaged magnetic stress (m xy ) 
roughly equals the stress from the horizontally-averaged fields, — (B x ) (B y ) /4.7T, except 
for two periods of a few orbits around 87 and 117 years. 
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The horizontally-averaged stress varies around zero, while in the ideal-MHD case it is 
strictly positive (figure [8]). 

The mean ratio of magnetic to hydrodynamic stresses is just 1.9, intermediate between 
the ratios in the dead zone and active layers during the early part of the calculation. 

The mean midplane pressure in the vertical component of the magnetic field is 125 times 
smaller than that in the azimuthal component, contrasting with the r atios around 30 



found in run Ml and also typical of past MRI turbulence calculations (jMiller fc Stone 



20001 ) . The relative strength of the azimuthal field suggests a generation mechanism 



distinct from the MRI. 

The mean azimuthal magnetic field evolves mostly through shear acting on the mean 
radial field. If shear were the only effect, the azimuthal field at a time t > t would be 



Bf{t) = By(t ) - [ jWdt 

JtQ z 



(11) 



The mean azimuthal field is plotted together with By in figure [9j The departure from 
the shear curve is shown over 20-year intervals by performing the integration from 
to = up to 10 years, from t = 10 up to 30 years, from to = 30 up to 50 years and 
so on. The field closely tracks the variation expected from shear throughout run VI, 
while the two curves are more weakly correlated in the ideal-MHD version Ml. 

These results indicate the magnetic fields evolve as follows. The fields delivered to the interior 
suffer Ohmic dissipation of small-scale structures, leaving a weak, almost uniform net radial 
component. Although the MRI is ineffective near the midplane due to the large resistivity, 
the shear produces a net azimuthal field, leading to an accretion stress that can be positive 
or negative. However the stress is negative for only short periods and when time-averaged 
is positive, driving accretion (figure [5]). The shear-generated fields are expelled both up 
and down, leading to correlated activity in the two surface layers unlike run Ml (figure [6]). 
The stresses and the azimuthal field strength peak about every te n orbits. Azimuthal field 



varia tions with similar periods occurred in MHD calculations by iBrandenburg fc Donner 



(119971 ) . The variations depend little on the domain height provided the box extends to 



l evels where the ma gnetic pressure exceeds the gas pressure, according to ideal MHD tests 



([Turner et al.l 120061 ). Occasionally when the midplane azimuthal field is strongest it erupts, 
leading to the production of a large midplane radial field as at 87 and 117 years, and then is 
destroyed within a few orbits (figure [9]) . The magnetic energy dissipated in these outbursts 
per unit area per unit time is comparable to the radiative flux from the surface of the viscous 
model disk with the same stress and mass column, suggesting such events could be directly 
observed. 
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Fig. 8. — Midplane magnetic accretion stresses at 1 AU in the ideal MHD calculation Ml 
(top panel) and the resistive MHD calculation VI (bottom panel). Stresses due to the entire 
magnetic field are shown by thick curves, those due to small-scale fields alone by thin red 
curves. 
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Fig. 9. — Midplane magnetic fields versus time at 1 AU in the ideal MHD calculation Ml (top 
panel) and the resistive MHD calculation VI (bottom panel). The azimuthal component is 
shown by thick solid curves, the radial component by dotted curves. Both are horizontally- 
averaged and the radial component is multiplied by ten so its variations can be read on 
the same scale. The azimuthal field that would result from shear acting on the radial field 
(eq. [HI) is indicated by thin red curves. 
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4.2. Turbulence at 5 AU 

Turbulence fills the entire disk thickness in the shearing-box calculations at 5 AU that 
have gas-phase magnesium abundances 10 -6 Solar and greater. No dead zone is present and 
the accretion stresses have almost identical mean vertical profiles in the cases M5, F52 and 
F56 (figure [TU1) . By contrast, the run F58 with magnesium abundance 10~ 8 Solar has a dead 
zone extending two scale heights from the midplane. The midplane magnetic dissipation 
time per scale height in F58 is just 34 orbits, compared with 630 orbits in run F56 (figure H]). 
Recombination at the midplane takes only three orbits and is much faster than even the 
mixing in the ideal MHD run (figure [3D , so turbulent mixing will not expunge the dead zone 
in this case. 



5. ACTIVITY CRITERIA 



Simple criteria for the presence of turbulence driven by the MRI are extremely useful 
in developing protostellar disk models without the investment of time needed for three- 
dimensional MHD calculations. Two such criteria are tested in this section. 

An Ohmic resistivity quenches the linear growth of the MRI if the magnetic fields diffuse 
across the characteristic wavelength A c faster than the perturbations grow. The diffusion time 
is the wavelength squared divided by the resistivity, yielding a criterion for instability on 
vertical fields 



Lu = 



> 1 



(12) 



in agreement with linear analyses (jjinlll996l ; ISano &: Miyamalll999l ). The left-hand side of 
eq. [T2]is the dimensionless Lundquist number. In the fully-developed turbulence, the growth 
is most easily halted by diffusion perpendicular to the midplane because on average the field 
reverses direction most often per unit length along the z-direction. Num erical unstratified 
shearing-box results indica te that eq. [12] applies to the turbulence also (ISano fc Inutsuka 
200ll ; ISano fc Stonell2002bl ). 



The second activity criterion is constructed using the diffusion time across the density 
scale height H. The resulting magnetic Reynolds number Re = c^/(rjQ) is independent of 
the field strength and so the value required for instability depends on the fields present. 
Turbulence occurred only at Reynolds numbers greater than 10 4 with zero net magnetic flux 
and at Reynolds numbers a s low as 100 w i th a net vertical magnetic flux, in unstratified 
shearing-box calculations by iFleming et all (120001 ) . 



The edges of the dead zone according to the two activity criteria are shown in figure [6j 
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Fig. 10. — Mean accretion stress versus height in four MHD calculations centered at 5 AU. 
The magnetic stress is shown by solid lines, the hydrodynamic stress by dotted lines. The 
four panels show results from an ideal MHD calculation (top left) and calculations with a 
fixed resistivity profile and gas-phase magnesium abundances 0.01 (top right), 10~ 6 (bottom 
left) and 10~ 8 times Solar (bottom right). The results are averaged from 10 to 100 orbits. 
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In run VI with the reaction network, the Reynolds number is less than 10 4 throughout the 
domain after about 60 years and greater than 100 everywhere from 85 years on. By the zero 
net flux criterion Re > 10 4 the dead zone fills the domain, which is clearly not the case. 
On the other hand by the vertical field criterion Re > 100, the dead zone is completely 
absent and MRI activity extends to the midplane, contradicting the evidence outlined in 
section 14.11 It may be that neither of these criteria applies to our calculations, which have 
zero net vertical but non-zero net radial and azimuthal magnetic fluxes. An intermediate 
critical Reynolds number of 1000, corresponding to a dissipation time of 160 orbits per scale 
height, more nearly approximates the boundaries of the MRI active regions, but fails to 
respond to the changes in magnetic field strength over the course of the calculation. In 
contrast to these three possible critical Reynolds numbers, as shown in the bottom panel of 
figure M a Lundquist number unity accurately traces the edges of the MRI turbulent layers 
at all times in run VI. 

The activity criteria are similar at 5 AU. The mean midplane Reynolds number is 
6 x 10 4 in run F52, 4000 in run F56 and 200 in run F58. Only the last of these has a dead 
zone (figure [TUT) . consistent with a critical Reynolds number of 1000 for this situation. The 
Lundquist number is greater than unity at all heights in the first two cases and greater than 
unity in run F58 only outside the MRI dead zone, again obeying eq. [121 

Overall, the Reynolds number suffers from the difficulty that its critical value must 
be determined for each magnetic field configuration using 3-D MHD calculations, while the 
Lundquist number always has a critical value unity. The Lundquist criterion is preferable 
because it explicitly shows the magnetic field dependence. Neither dimensionless number 
bypasses the basic uncertainty regarding the net magnetic flux, which is an input parameter 
for the shearing box calculations. However once a disk model and field are chosen, the 
Lundquist criterion eq. [T21 immediately gives the size of the dead zone. 



6. CONCLUSIONS 



We investigated whether the dead zone is eliminated from the minimum-mass Solar 
nebula under favorable conditions. Cosmic rays are assumed to reach the disk surface, 
recombination on grains is absent and metals are found in the gas phase. The most restrictive 
of these assumptions is the lack of small grains. Particles 0.1 fim in radius can significantly 
affect the thickness of the dead zone at 1 AU ey en with dust-to-gas mass ratios 100 times 
lower than the interstellar value (jSano et al.ll2000l ). The gas-phase metal abundance is a less 
important parameter and can be several orders of magnitude below its nominal value with 
little change in the dead zone size. 
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Timescale estimates and numerical resistive MHD calculations both show that the whole 
of the disk thickness at 5 AU is active if the gas-phase magnesium abundanc e is greater than 



10 6 ti mes Solar, consistent with results at a similar column density from lllgner fc Nelson 



(l2006al ). The problem is more severe at 1 AU, where few ionizing particles penetrate to 
the midplane if the column density is 1700 g cm -2 . In chemical equilibrium the dead zone 
thickness declines with increasing gas-phase magnesium abundance but asymptotes at values 
above 10 -4 times Solar with most of the mass column still magnetically inactive. However 
ionized gas can reach the interior through turbulent mixing, because the recombination and 
mixing timescales are similar. In an MHD calculation including time- dependent ionization 
chemistry, free electrons arrive at the midplane fast enough to maintain a weak coupling to 
the magnetic fields. The midplane remains linearly stable against the MRI because of the 
large resistivity, but receives magnetic fields generated in the surface layers. The midplane 
fields are smooth while those in the surface layers are tangled. The boundary between the 
two lies at the height where the MRI is marginally stable and the Lundquist number is 
unity. Radial fields of one sign preferentially escape from the disk surfaces, leaving the 
midplane with a net radial magnetic flux of the opposite sign and leading to the generation 
of azimuthal fields through shear. The resulting accretion stress varies around zero, but on 
average is positive and just a few times less than the stress in the active layers. 

These results demonstrate that the dead zone can be eradicated and the midplane gas 
flow toward the star even if the MRI occurs only in the surface layers. Magnetic forces 
can provide angular momentum transfer throughout the planet formation region once small 
grains have been removed. However the midplane flow at 1 AU is fundamentally differ- 
ent from the usual picture. It is magnetorotationally stable and has magnetic stresses that 
are uniform on scales up to the disk thickness. The presence of internal stresses yields 
a greater overall accretion rate than the basic layered model, permitting a stronger de- 
pendence on the stella r mass as required by brown dwarf and T Tauri star measurements 



( iHartmann et al.ll2006l ). Time variability in our calculations is due to fluctuations in the 



magnetic field strength. Longer-period accretion variability can arise f rom the pileup of ma - 



terial reaching the dead zone, even in the presence of midplane stresses (jWunsch et al.ll2006l ); 
these effects are not included in our calculations. The smooth internal magnetic fields have 
pressures up to 5% of the gas pressure and are strong enough to modify the orbit al migra- 



tion of Earth-mass protoplanets through their effects on the corotation resonance (ITerquem 



20031 ). Other potentially important effects in the dead zone are the transport of magnetic 
flux over large radial distances and the interactions between disk annuli carrying opposing 
magnetic fields. These can be explored with future global MHD calculations. 



Portions of this work were carried out at the Jet Propulsion Laboratory, California 
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Institute of Technology using the JPL Supercomputing and Visualization Facility. Support 
was provided by the JPL Research and Technology Development and NASA Solar System 
Origins Programs. 

REFERENCES 

Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in "Cosmic Abundances as Records of 
Stellar Evolution and Nucleosynthesis", ASP Conf. Ser., vol. 336, ed. T. G. Barnes 
III & F. N. Bash (San Francisco: Astronomical Society of the Pacific), 25 

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 

Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 

Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235 

Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 

Brandenburg, A., k Donner, K. J. 1997, MNRAS, 288, L29 

Casein, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 

Fitzpatrick, E. L. 1997, ApJ, 482, L199 

Fleming, T., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464 

Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 

Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751 

Fromang, S., Terquem, C, & Balbus, S. A. 2002, MNRAS, 329, 18 

Gammie, C. F. 1996, ApJ, 457, 355 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 
Hawley, J. F., & Stone, J. M. 1998, ApJ, 501, 758 

Hartmann, L., D'Alessio, P., Calvet, N., & Muzerolle, J. 2006, ApJ, 648, 484 

Hayashi, C, Nakazawa, K., & Nakagawa, Y. 1985, in: Black, D. C, & Matthews, M. S., 
Editors. Protostars and Planets II (Univ. of Arizona Press: Tucson), 1100 

Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 



-24- 



Ilgner, M., & Nelson, R. P. 2006a, A&A, 445, 205 
Ilgner, M., & Nelson, R. P. 2006b, A&A, 445, 223 
Inutsuka, S., & Sano, T. 2005, ApJ, 628, L155 
Jin, L. 1996, ApJ, 457, 798 

Matsumura, S., & Pudritz, R. E. 2006, MNRAS, 365, 572 
Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398 
Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 
Oppenheimer, M., & Dalgarno, A. 1974, ApJ, 192, 29 
Pneuman, G. W., & Mitchell, T. P. 1965, Icarus, 4, 494 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, "Numerical 
Recipes in C" , 2nd ed (Cambridge: Cambridge University Press) 

Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93 

Salmeron, R., & Wardle, M. 2005, MNRAS, 361, 45 

Sano, T., & Inutsuka, S. 2001, ApJ, 561, L179 

Sano, T., & Miyama, S. M. 1999, ApJ, 515, 776 

Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 

Sano, T., & Stone, J. M. 2002a, ApJ, 570, 314 

Sano, T., & Stone, J. M. 2002b, ApJ, 577, 534 

Shakura, N. L, & Sunyaev, R. A. 1973, A&A, 24, 337 

Sofia, U. J., Cardelli, J. A., & Savage, B. D. 1994, ApJ, 430, 650 

Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753 

Stone, J. M., & Norman, M. L. 1992b, ApJS, 80, 791 

Terquem, C. E. J. M. L. J. 2003, MNRAS, 341, 1157 

Turner, N. J., Willacy, K., Bryden, C, & Yorke, H. W. 2006, ApJ, 639, 1218 



-25 - 



Umebayashi, T. 1983, Progr. Theor. Phys., 69, 480 

Wardle, M. 1999, MNRAS, 307, 849 

Wardle, M., & Ng, C. 1999, MNRAS, 303, 239 

Weidenschilling, S. J., & Cuzzi, J. N. 1993, in: Levy, E. H., & Limine, J. I., Editors. 
Protostars and planets III (Univ. of Arizona Press: Tucson), 1031 

Wiinsch, R., Gawryszczak, A., Klahr, H., & Rozyczka, M. 2006, MNRAS, 367, 773 



This preprint was prepared with the AAS IATgX macros v5.2. 



